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ABSTRACT 

The  optimization  method  proposed  in  this  paper  is  for  determining  open  boundary  conditions  from  interior 
observations.  Unknown  open  boundary  conditions  are  represented  by  an  open  boundary  parameter  vector  (B), 
while  known  interior  observational  values  are  used  to  form  an  observation  vector  (O).  For  a  hypothetical  B* 
(generally  taken  as  the  zero  vector  for  the  first  time  step  and  as  the  optimally  determined  B  at  the  previous  time 
step  afterward),  the  numerical  ocean  model  is  integrated  to  obtain  solutions  (S*)  at  interior  observation  points. 

The  root-mean-square  difference  between  S*  and  O  might  not  be  minimal.  The  authors  change  B*  with  different 
increments  SB.  Optimization  is  used  to  get  the  best  B  by  minimizing  the  error  between  O  and  S. 

The  proposed  optimization  method  can  be  easily  incorporated  into  any  ocean  models,  whether  linear  or 
nonlinear,  reversible  or  irreversible,  etc.  Applying  this  method  to  a  primitive  equation  model  with  turbulent 
mixing  processes  such  as  the  Princeton  Ocean  Model  (POM),  an  important  procedure  is  to  smooth  the  open 
boundary  parameter  vector.  If  smoothing  is  not  used,  POM  can  only  be  integrated  within  a  finite  period  (45 
days  in  this  case).  If  smoothing  is  used,  the  model  is  computationally  stable.  Furthermore,  this  optimization 
method  performed  well  when  random  noise  was  added  to  the  “observational”  points.  This  indicates  that  real¬ 
time  data  can  be  used  to  inverse  the  unknown  open  boundary  values. 


1.  Introduction 

One  of  the  difficult  problems  in  shallow-water  mod¬ 
eling  is  the  uncertainty  of  the  open  boundary  condition 
(OBC).  At  open  boundaries  where  the  numerical  grid 
ends,  the  fluid  motion  should  be  unrestricted.  Ideal  open 
boundaries  are  transparent  to  motions.  Two  approaches, 
local  type  and  inverse  type,  are  available  for  determining 
the  OBC.  The  local-type  approach  determines  the  OBC 
from  the  solution  of  the  governing  equations  near  the 
boundary.  The  problem  now  becomes  selection  from  a 
set  of  ad  hoc  OBCs.  Since  any  ad  hoc  OBC  will  intro¬ 
duce  inaccuracies  into  a  numerical  solution  (Chapman 
1985),  it  is  important  to  choose  a  best  one  from  ad  hoc 
OBCs  for  a  particular  ocean  model.  Using  a  barotropic 
coastal  ocean  model.  Chapman  (1985)  evaluated  several 
of  the  most  used  ad  hoc  OBCs  (clamped,  sponge,  ra¬ 
diation)  and  found  that  the  best  OBC  consists  of  a 
sponge  at  the  outer  edge  of  the  model  domain  with  an 
Orlanski  radiation  condition  (Orlanski  1976)  and  that 
the  clamped  OBC  is  probably  the  worst  choice.  Apply¬ 
ing  these  results  to  other  ocean  models  needs  further 
investigation.  The  local  approach  suffers  drawbacks  that 
may  restrict  its  use:  no  observational  data  considered 
and  the  ill-posedness  of  the  primitive  equations  model 
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with  ad  hoc  OBC,  that  is,  it  is  hard  to  prove  the  existence 
of  a  unique  solution  (Bennett  1992;  Oliger  and  Sund- 
strom  1978).  To  improve  the  local  approach  by  using 
observations  at  open  boundaries,  Shulman  and  Lewis 
(1995)  proposed  a  method  for  determining  OBCs  of  the 
shallow-water  model.  Their  method  is  based  on  the  in¬ 
tegration  of  governing  equations  forward  in  time  and 
the  selection  of  OBCs  via  a  specific  inverse  problem 
that  minimizes  a  measure  of  difference  (energy  flux) 
between  the  values  of  observed  and  predicted  variables 
at  open  boundaries.  Thus,  their  method  helps  us  to  select 
the  proper  ad  hoc  OBC  by  using  observations  at  the 
open  boundaries. 

Without  any  ad  hoc  OBCs,  the  inverse-type  approach 
can  determine  the  OBC  from  a  “best”  fit  between  model 
solutions  and  interior  observations.  The  most  popular 
scheme  for  this  approach  is  an  adjoint  method,  which 
consists  of  four  elements:  set  of  control  parameters  or 
control  vector  (e.g.,  the  unknown  OBC),  numerical 
ocean  model,  cost  function,  and  adjoint  equation.  The 
cost  function  is  usually  defined  by  the  difference  be¬ 
tween  observations  and  their  model  counterparts.  The 
adjoint  equation  is  derived  from  minimizing  the  cost 
function  with  respect  to  the  control  parameters.  Using 
an  adjoint  method,  the  initial-value  problem  of  ocean 
model  with  the  OBC  becomes  integration  of  both  the 
governing  equations  and  the  equation  for  the  control 
parameters  forward  and  backward  in  time.  For  a  com¬ 
prehensive  discussion  of  the  adjoint  method,  the  reader 
is  referred  to  the  numerous  papers  on  that  subject,  for 
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example,  Seiler  (1993).  The  advantage  of  using  the  ad¬ 
joint  method  is  the  well  posedness  and  the  use  of  ob¬ 
servational  data.  Seiler  (1993)  successfully  determines 
the  unknown  OBCs  for  a  quasigeostrophic  ocean  by 
using  the  adjoint  method.  The  disadvantages  that  may 
restrict  its  use  are  the  requirement  of  large  amounts  of 
computer  time  and  memory;  problems  of  stable  inte¬ 
gration  of  the  adjoint  equation;  ocean-model  dependen¬ 
cy  of  the  adjoint  equation;  and  difficulty  in  deriving  the 
adjoint  equation  when  the  model  contains  rapidly  chang¬ 
ing  processes,  such  as  ocean  mixed  layer  dynamics. 

We  propose  a  simplified  method  that  overcomes  the 
disadvantage  of  the  current  inverse-type  approach.  This 
method  can  determine  OBCs  of  any  ocean  model  (i.e., 
a  universal  method)  from  interior  observations.  The  es¬ 
sence  of  the  method  is  to  seek  the  relationship  among 
three  vectors:  open  boundary  parameter  vector  B,  ob¬ 
servation  vector  O,  and  solution  vector  S.  If  B  is  given, 
we  can  integrate  the  numerical  ocean  model  and  get  the 
solution  vector  S.  If  B  is  unknown,  the  optimization 
method  is  used  to  determine  B  by  minimizing  the  root- 
mean-square  (rms)  difference  between  O  and  S. 

2.  Optimization  method 

a.  Three  vectors 

If  r  denotes  the  position  of  any  point  along  the  open 
boundary  T,  the  boundary  values  of  any  variable  17  is 
a  function  of  r,  yfb)  =  rfy,(r).  Let  f\(r),  f2(r),  .  .  .  ,f„(r) 
be  a  series  of  known  basis  functions.  We  expand  the 
function  rfb)(r)  into 

n 

Vlh,(r)  =  2  hifi(r)-  (1) 

1=1 

Thus,  the  determination  of  the  open  boundary  condition 
rfh,(r)  becomes  the  determination  of  a  set  of  parameters 
bx,  b2,  ■  ■  ■  ,  bn.  The  ^-dimensional  vector. 


is  called  the  boundary  parameter  vector. 

Assume  that  there  are  m  observations  (Ou  02,  .  .  . , 
Om )  located  at  the  interior  (Fig.  1).  An  m-dimensional 
vector  O  can  be  constructed  by 
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Fig.  1.  Can  open  boundary  condition  be  determined  by  interior 
values? 


Fig.  2.  Flowchart  showing  the  optimization  method  for  determin¬ 
ing  open  boundary  vector  B  for  the  kth  time  step. 
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Time  Step 

Fig.  3.  Smooth  of  the  open  boundary  correction  vector  by  the  linear 
regression  method.  Circles  indicate  smoothed  data  and  asterisks  de¬ 
note  unsmoothed  data. 

which  is  called  the  observation  vector.  Notice  that  the 
observational  variable  is  not  necessarily  the  same  di¬ 
mension  as  the  variable  at  the  open  boundary.  If  B  is 
given,  we  can  solve  the  dynamic  system  and  obtain  the 
solution  S.  At  the  same  locations  where  the  observations 
take  place,  the  solutions  are  Slt  S2,  .  .  . ,  Sm,  which  form 
a  solution  vector 


70  45  90  135  180 

CROSS  SHORE  (km) 


x 

Fig.  4.  Integration  domain  and  lateral  boundary  conditions. 


b.  Optimization  method  for  determining  B 

Ocean  model  performance  can  be  measured  by  the 
rms  error 

“|  1/2 

j  m 

I  =  -  2  (Sj  -  oy  .  (5) 

m  ,  1 

The  vector  S  depends  on  B.  Change  of  the  boundary 
vector  B  (boundary  conditions)  leads  to  a  change  of  S 
(solutions).  Inversely,  we  may  determine  B  by  mini¬ 
mizing  /; 

31 

—  =  0,  i  =  1,  2,  .  . . ,  n.  (6) 

db, 

Substitution  of  (5)  into  (6)  leads  to  a  set  of  n  equations 
implicitly  solvable  for  bu  b2,  ... ,  bn, 

m 

2  CS,  -  0J)R,J  =  0,  *  =  1,  2 . n,  (7) 

j=  1 

where 

dS: 

Rii  =  — i  =  1,  2,  . . . ,  n;  j  =  1,  2,  .  . . ,  m 

,j  db,  j 

(8) 

are  components  of  a  n  X  m  Jacobian  matrix  R  =  { Rlf ) . 


From  a  first-guess  boundary  vector  B*,  a  solution 
vector  S*  is  obtained  by  solving  the  numerical  ocean 
model.  The  rms  between  S*  and  O  might  not  be  min¬ 
imal.  We  update  the  boundary  parameter  vector  com¬ 
ponents  by  increments  {<5/?,  |  i  =  1,  2, ... ,  n}.  and 
therefore  components  of  the  solution  vector  become 


fl  f2  13  f4  f5 


0  .5  10  .5  10  .5  1-10  10  1  2 

Basis  Functions 

Fig.  5.  Five  basis  functions  of  the  open  boundaries,  used  in  the 
Csanady’s  shelf  model. 


726 


JOURNAL  OF  ATMOSPHERIC  AND  OCEANIC  TECHNOLOGY 


Volume  14 


0  30  60  90  120  150  180 


0  30  60  90  120  150  180 


X-location  (km) 

Fig.  6.  Sensitivity  of  relative  errors  on  the  cross-coastal  location 
of  the  observation  points:  (a)  interior  error  and  (b)  open  boundary 
error. 


Sj  =  Sf  +  2  RijSbj  +  high-order  terms.  (9) 

i=i 

Substituting  (9)  into  (7)  and  neglecting  higher-order 
terms  leads  to  a  set  of  n  linear  algebraic  equations  for 
{»,}, 

n 

2  PuSb,  =  d„  i  =  1,2,...,  n,  (10) 
/=  1 

where 

m  m 

p„  -  2  RA,  d,  -  2  Ri]ioj  -  s*y, 

J= 1  (11) 

i  =  1,  2,  . .  . ,  n;  l  =  1,  2,  . . . ,  n. 

Both  0;  and  S*j  are  known  quantities.  Therefore,  the 
linear  algebraic  equations  (10)  have  definite  solutions 
when  the  Jacobian  matrix  { )  is  determined  and 

det{P,;}  7^  0.  (12) 

c.  Determination  of  the  Jacobian  matrix  by  a 
multiperturbation  method 

As  soon  as  the  Jacobian  matrix  is  obtained,  we  can 
solve  (10)  to  get  boundary  parameter  vector  corrections 


(a) 


01  1  •  - - — 

0  200  400  600  800  1000  1200  1400  1600  1800  2000 

x  (km) 


Fig.  7.  Two  areas  for  the  POM  integration:  (a)  domain  with  three 
rigid  and  one  open  boundaries  and  (b)  a  double  size  domain  with 
four  rigid  boundaries. 


Wind  Stress  (mA2/sA2)x  ^  q-4 


Fig.  8.  Pseudo-zonal  wind  stress  (m2  s  2)  used  in  the  POM  inte¬ 
gration. 
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(a)  (b) 
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Basis  Functions 
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Fig.  9.  Basis  functions  and  corresponding  amplitudes:  (a)  fi(y),  (b) 
f2(y),  and  (c)  bi(t)  and  b2(t). 
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0 
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0 

1 

and  perturb  the  first  guess  at  each  time  step  by  a  small 
fraction  of  these  unit  vectors  c®: 

B®  =  B*  +  efcc®,  i  =  1,  2,  .  .  . ,  n,  (16) 

where  eb  is  a  small  positive  number.  Using  the  n  +  1 
sets  of  boundary  parameter  vectors  B*,  B(1),  B(2),  .  .  .  , 
B‘"),  we  obtain  n  +  1  sets  of  solution  vectors  for  that 
time  step,  S*,  S(1),  S,2),  .  .  .  ,  S<").  The  Jacobian  matrix 
(8)  denotes  the  rate  of  change  of  the  solution  vector 
with  the  boundary  parameter  vector;  therefore,  it  can 
be  computed  by 

1 

R  =  —  [e«e<2>  •  •  ■  e(n)],  (17) 

where  the  m-dimensional  vectors  eu),  et2),  .  .  .  ,  e1”’  are 
defined  by 

e®  =  S®  -  S*.  *  =  1,2,...,  n.  (18) 


d.  Iteration  process 

After  the  algebraic  equation  (10)  is  solved,  the  bound¬ 
ary  parameter  correction  vector  SB  =  ( 8b , ,  8b2,  .  .  .  , 
8b„)  is  obtained.  We  replace  B*  by  B*  +  SB  and  repeat 
the  process  (Fig.  2)  until  reaching  a  certain  criterion 

|B*|  ’ 

where 


8bi  (i  =  1,2,...,  it).  Therefore,  the  determination  of 
the  Jacobian  matrix  is  the  key  issue.  We  propose  a  sim¬ 
ple  multiperturbation  method  on  the  first-guess  bound¬ 
ary  parameter  vector  B*  to  get  the  Jacobian  matrix.  The 
vector  B*  at  the  initial  state  ( 1  =  0)  of  any  numerical 
ocean  model  is  set  to  be  a  zero  vector. 


0 

0 


B*(0)  = 


0 


(13) 


and  the  first  guess  for  B  at  time  step  k  is  set  to  be  the 
same  as  the  determined  boundary  parameter  vector  at 
time  step  k  —  1, 

B*(£)  =  B (k  -  1),  k  =  1,  2, _  (14) 

We  define  a  set  of  H-dimensional  unit  vectors 


I B*  |  -  l-^br 


I  SB  | 


(19) 


and  e  is  a  small  positive  number  (user  input).  As  soon 
as  the  inequality  (19)  is  satisfied,  the  iteration  stops  and 
the  final  B*  becomes  the  optimal  boundary  parameter 
vector  for  the  next  time  step. 


e.  Reference  model 

A  model  with  a  given  boundary  condition  (called 
reference  boundary  condition)  is  needed  for  error  es¬ 
timation.  We  run  the  model  with  the  given  reference 
boundary  condition  and  obtain  the  solution  for  interior 
points,  O  =  (O „  02,  .  .  .  ,  (9,„),  which  are  taken  as  “ob¬ 
servations.”  We  also  expand  the  reference  boundary 
values  into  basis  functions  (1)  to  obtain  the  reference 
boundary  vector  B(ref)  =  (b'f,  b2{,  .  .  .  ,  bf).  Such  a  mod¬ 
el  (with  known  boundary  conditions)  is  called  a  refer¬ 
ence  model. 
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Fig.  10.  Time  series  of  open  boundary  values  of  the  domain  A.  The  dots  represent  the 
values  computed  from  the  basis  functions. 


f.  Smoothing  8 B 

Both  observational  (instrument)  and  computational 
(numerical)  errors  perturb  the  values  of  B.  When  the 
reference  model  results  are  used,  there  is  no  instrument 
error  in  observations.  But  the  computational  errors  still 
cause  errors  in  B.  For  prognostic  ocean  models,  the 
errors  of  B  feed  back  to  the  next  step  model  integration. 
Error  accumulation  in  each  time  step  may  cause  com¬ 
putational  instability.  We  use  linear  regression  for 
smoothing  SB  to  reduce  high-frequency  modes.  The 
procedure  is  as  follows:  1)  For  the  first  11  time  steps, 
we  do  not  do  any  smoothing,  and  after  the  11th  time 
step,  we  use  the  smoothing  technique.  2)  For  each  time 
step  k  (k  >  11),  we  maintain  smoothed  values  of  each 
component  at  the  previous  11  steps,  8bf\k  —  11), 
8b)"(k  —  10),  .  .  .  ,  8b)s){k  —  1).  Here  the  superscript  s 
indicates  smoothed  value.  3)  We  fit  a  linear  regression 
to  fit  the  12  values:  8b)s)(k  -  11),  8b<f(k  -  10),  .  .  . , 
8b)s)(k  —  1),  and  8bj(k).  4)  We  obtain  the  value  at  the 
time  step  k  from  the  regression,  Bb'/Hk).  5)  We  replace 
8bj(k)  by  8b(f(k).  Such  a  treatment  (Fig.  3)  will  filter 
out  high-frequency  noise  in  the  derived  open  boundary 
parameter  vector. 


F(80)  = 


V27 tc 


exp 


(soy 


2a2 


(20) 


where  80  is  a  random  variable  with  a  zero  mean  and  a 
standard  deviation  of  cr. 


h.  Relative  errors 


For  a  given  observation  vector  O  and  unknown  open 
boundary  condition,  we  use  the  optimization  method  to 
obtain  the  values  at  the  open  boundary,  B  =  (/?,,  b2, 
. .  .  ,  bn).  Then  we  integrate  the  same  model  with  the 
computed  open  boundary  condition  and  get  the  solutions 
at  the  observational  points,  S  =  (5),  S2,  .  .  . ,  Sm).  The 
relative  errors. 


Ew 


2  W  ~  b\ 

I \b?\ 


Em 


1\Oj-Sj\ 


(21) 


measure  the  validity  of  the  optimization  method.  The 
smaller  the  E(m  and  Eto>,  the  better  the  optimization 
method. 


g.  Random  noises  added  on  “observations” 

When  observations  are  taken  from  the  reference  mod¬ 
el  solution,  there  is  no  observational  (instrument)  error. 
In  order  to  see  the  effects  of  observational  error  on  the 
determination  of  open  boundary  conditions,  we  add  a 
Gaussian-type  random  variable  to  each  observation  at 
each  time  step.  The  probability  distribution  function  is 
given  by 


3.  Linear  ocean  models 

a.  Jacobian  matrix 

For  any  linear  ocean  model,  the  relationship  between 
{Sj}  and  {Sbj}  becomes  linear.  There  are  no  high  order 
terms  in  (9).  The  Jacobian  matrix  can  be  easily  obtained 
by  letting 
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Fig.  1 1 .  The  velocity  fields  from  the  reference  model  results  after 
(a)  30  day's  run  and  (b)  90  day’s  run. 


E 


Fig.  12.  Randomly  picked  observational  points  for  the  POM  model. 


b.  Example — Csanady’s  shelf  model 

We  use  a  steady-state  shelf  circulation  model  pro¬ 
posed  by  Csanady  (1978)  as  an  example  to  show  how 
to  obtain  the  Jacobian  matrix  R  and  to  determine  the 
open  boundary  conditions.  Consider  a  long  and  straight 
coastline  with  coordinates  such  that  the  y  axis  coincides 
with  the  coast  and  positive  x  points  offshore.  The  water 
depth  is  a  function  only  of  offshore  distance,  that  is,  h 
=  h(x).  As  pointed  out  by  Csanady  (1978),  this  model 
is  appropriate  for  simulating  mean  flow  along  such  a 
coast  driven  by  mean  wind  stress. 

Under  the  assumptions  of  homogeneous  water  and 
linear  bottom  friction  (i.e.,  proportional  to  depth-aver¬ 
aged  velocity),  the  dynamical  equations  for  depth-av¬ 
eraged  velocities  ( u ,  v)  without  wind  forcing  become 

dri 

~fv  =  -g-r-  (22) 

OX 


3  ri  rv 

f“  *  _ss7  ”  T, 


B(uh)  d(vh) 

dx  dy 


(23) 

(24) 


for  any  time  step.  The  algebraic  equation  (9)  becomes 


where  17  is  the  surface  elevation,  r  the  bottom  resistance 
coefficient,  /  the  Coriolis  parameter,  and  g  the  gravity. 
Eliminating  u  and  v  from  (22)-(24)  leads  to  a  single 
equation  for  the  surface  elevation  17: 


d2V  +  fdhdp  =  Q 
dx2  r  dx  dy 


(25) 


Wang  (1982)  applied  the  Csanady  model  to  a  region 
(Fig.  4)  with  a  combined  flat  shelf  and  steep  slope, 
depicted  by 


h(x)  = 


0  <  x  <  x0 
x0  =£  x  <  xu 
X  >  xu 


RIV‘>  =  e‘" 


i  =  1,2 ,...,«. 


0.02  +  10-3x, 

0.16  +  0.050r  —  x0), 

2.0, 
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(a) 


(c) 


Fig.  13.  Components  of  the  open  boundary  parameter  vector  ob¬ 
tained  by  the  optimization  method  for  (a)  case  1.  (b)  case  2,  and  (c) 
case  3. 


where  x0  =  140  km  and  x,  =  180  km  are  the  offshore 
location  of  the  shelf  break  and  slope  edge,  respectively. 
Following  Wang  (1982),  we  chose /=  10~4  s~\  and  r 
=  0.1  cm  s  ',  and  use  the  boundary  condition  at  the 
coast  (x  =  0), 


dr]  r  dr] 

dy  fh  (x0)  dx 


=  0. 


(26) 


We  consider  the  pure  open  ocean  forcing  case;  that  is, 
there  is  no  inflow  at  the  backward  boundary  (v  =  0), 


77  =  0,  at  y  =  0.  (27) 

Equation  (25)  has  the  form  of  the  one-dimensional 
heat-conduction  equation,  with  negative  y  playing  the 
role  of  time.  We  solve  (25)  numerically  by  the  Gaussian 


elimination  method  under  the  boundary  conditions 
(26)-(27)  to  obtain  the  first-guess  solution  at  the  interior. 
The  spatial  increments  are  Ax  =  2  km,  A y  =  10  km. 
We  use  the  five  basis  functions  fx(y),  /2(y),  ■  ■  ■ ,  fs(y ) 
(Fig.  5)  for  7]a,).  The  first  guess  for  B  is  the  zero  vector. 
For  simplicity  and  no  loss  of  generality,  we  choose  five 
equally  spaced  points  on  an  interior  line  paralleling  to 
the  y  axis  as  observation  points.  In  this  case. 


1  0  0  0  0 
0  10  0  0 
0  0  10  0 
0  0  0  1  0 
0  0  0  0  1 


(28) 


is  a  unit  matrix.  The  Jacobian  matrix  then  becomes 

R  =  [e(ne<2,e(3,e(4'e(5)].  (29) 


c.  Error  estimation 

For  simplicity,  we  choose  Wang’s  (1982)  solution 
along  with  the  open  boundary  condition  as  a  reference 
model  for  the  error  estimation.  Five  observation  points 
are  equally  spaced  and  located  along  lines  paralleling 
to  the  y  axis.  These  lines  are  represented  by  their  x 
coordinate.  Wang’s  (1982)  solution  at  the  five  points  is 
taken  as  observations,  O  =  (O,,  02,  03,  04,  05),  and 
the  Wang’s  open  boundary  vector  is  denoted  by  B(ref)  = 
(b'f,  b'f,  b'f,  b'f,  (?5et).  Since  the  model  is  steady  state, 
we  did  not  smooth  SB. 

The  five  observation  points  are  given  as  O  with  un¬ 
known  open  boundary  condition.  Using  the  optimization 
method,  we  obtain  values  at  the  open  boundary,  B  = 
(blt  b2,  b},  b4,  lx).  Then,  we  integrate  the  Csanady  model 
(24)  with  the  computed  open  boundary  condition,  and 
get  the  solutions  at  the  same  five  points,  S  =  (5),  S2, 
S3,  S4,  S5).  Since  the  observation  points  are  on  the  lines 
paralleling  the  y  axis  in  this  study,  E(B)  and  E'°>  depend 
only  on  x.  Surprisingly,  both  E(B)  and  E'°'  are  extremely 
small  (Fig.  6).  When  the  observation  points  are  chosen 
near  the  coast,  the  relative  errors  are  on  the  order  of 
10  6— 1 0  7.  At  the  shelf  break  (x  =  x0,  x0  =  140  km), 
EtB)  is  on  the  order  of  10~6.  Passing  the  shelf  break,  E(B) 
and  Eta)  decrease  very  fast  offshore  and  are  on  the  order 
of  10“ 16  near  the  open  boundary. 

4.  Nonlinear  ocean  model 

a.  Example — Princeton  Ocean  Model  (POM) 

We  apply  the  optimization  method  to  determine  the 
open  boundary  conditions  of  a  flat  bay  centered  at  35°N 
and  bounded  by  three  rigid  boundaries.  This  bay  ex¬ 
pands  1000  km  in  both  the  north-south  and  east-west 
directions.  The  northern,  southern,  and  western  bound¬ 
aries  are  rigid,  and  the  eastern  boundary  is  open.  The 
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x  (km) 


Cartesian  coordinate  system  is  chosen  with  the  origin 
at  the  southwest  corner.  The  x  axis  points  toward  the 
east,  and  the  y  axis  toward  the  north  (Fig.  7a).  The 
circulation  in  the  bay  is  modeled  with  the  Princeton 
Ocean  Model  (POM)  developed  by  Blumberg  and  Mel- 
lor  (1987).  POM  is  a  primitive  equation  model  with  a 
free  surface  and  a  level-2  turbulence  closure  scheme 
(Mellor  and  Yamada  1982).  A  description  of  the  model 
code  can  be  found  in  Mellor  (1991).  We  use  the  2D 
version  of  POM  to  illustrate  the  usefulness  of  the  op¬ 
timization  method  for  determining  the  open  boundary 
conditions. 

The  area  depicted  in  Fig.  7a  is  called  domain  A,  where 
the  boundary  conditions  are  known  at  the  three  rigid 
boundaries  (northern,  southern,  and  western),  and  un¬ 
known  at  the  eastern  boundary.  The  eastern  boundary 
of  domain  A  is  connected  to  a  mirror  image  of  domain 
A  (about  x  =  1000  km)  forming  a  closed  rectangular 
domain  (Fig.  7b),  called  domain  B.  The  POM  model 
was  integrated  from  the  following  initial  conditions. 


(b) 


x  (km) 


Fig.  14.  Horizontal  velocity  fields  after  the  30-day  run: 
(a)  case  1,  (b)  case  2,  and  (c)  case  3. 


(u,  V,  w)  =  0  (m  s-1),  T=  283(1  +  erfloo°)  (K), 

S  =  35  (psu),  (30) 

under  no  surface  heat  or  salinity  fluxes  and  zonal  surface 
pseudo-wind  stress  varying  with  latitude  (Fig.  8): 

—  =-10-4cos—  (m2  s-2).  (31) 

Po  Ly 

The  time  step  was  chosen  as  2  min.  The  horizontal 
resolution  was  50  km.  Bottom  stress  is  parameterized 
by  the  quadratic  drag  relation.  Horizontal  kinematic  vis¬ 
cosity  is  set  to  be  500  m2  s_1. 

b.  Boundary  vector 

We  integrate  POM  over  domain  B  with  four  rigid 
boundaries  (known  boundary  conditions)  from  the  ini¬ 
tial  conditions  (30)  and  surface  forcing  (31)  and  take 
the  solution  along  the  middle  of  domain  B  (x  =  1000 
km)  as  the  reference  open  boundary  condition  for  the 
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Fig.  16.  The  relative  boundary  vector  error  for  (a)  case  1,  (b)  case 
2,  and  (c)  case  3. 


Fig.  15.  Horizontal  velocity  fields  after  the  90  day’s  run:  (a)  case  2 
and  (b)  case  3. 


B(0  =  [*,(0,  b2m  (33) 


domain  A  integration.  The  velocity  at  x  =  1000  km  for 
the  domain  B  run  is  nearly  zonal, 

|v|  <C  \u\  at  x  =  1000  km. 

For  simplicity  and  no  loss  of  generality,  we  assume  zero 
latitudinal  velocity  at  the  eastern  open  boundary  for  the 
domain  A  integration.  Therefore,  our  task  here  is  to 
determine  the  zonal  velocity  at  the  eastern  open  bound¬ 
ary.  Following  the  procedure  depicted  in  section  2,  we 
expand  the  boundary  values  by  two  basis  functions/'^}’) 
and  f2(y)  (Figs.  9a,b), 

u^(y,  t )  =  bft)ffy)  +  b2(t)f2(y)  (32) 

with  two  temporally  varying  amplitudes  bft)  and  b2(t), 
as  shown  in  Fig.  9c.  After  plotting  the  time  series  (from 
day  5  to  day  30)  of  the  reference  open  boundary  values 
and  the  corresponding  expanded  values  (Fig.  10),  we 
find  that  the  boundary  parameters  can  be  well  repre¬ 
sented  by  the  two-dimensional  vector. 


c.  Reference  model 

We  choose  the  POM  model  solution  for  domain  A 
computed  with  the  reference  open  boundary  condition 
determined  by  the  two  time  series  bft),  b2(t )  (Fig.  9) 
as  a  reference  model.  We  ran  the  reference  model  for 
90  days  with  known  boundary  conditions  (three  rigid 
and  one  reference  open  boundary  condition).  Figure  1 1 
shows  the  horizontal  velocity  of  the  reference  model 
run  for  two  different  times  (30  day  and  90  day).  We 
will  use  the  reference  model  results  to  verify  the  opti¬ 
mization  method. 

d.  Three  cases 

A  random  number  generator  (Fortran  function,  Ranf) 
was  used  to  produce  random  disturbances  for  each  ob¬ 
servational  point  independently  with  mean  value  of  zero 
and  standard  deviation  of  0.01  m  s  '.  In  order  to  test 
the  performance  of  the  optimization  method,  we  ran 
three  cases:  1)  without  smoothing  on  SB  and  without 
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Fig.  17.  The  relative  error  at  the  observational  points  for  (a)  case 
1,  (b)  case  2,  and  (c)  case  3. 

random  noise  added  to  the  observations,  2)  with  smooth¬ 
ing  on  SB  and  without  random  noise  added  to  the  ob¬ 
servations,  and  3)  with  both  smoothing  on  SB  and  ran¬ 
dom  noise  added  to  the  observations. 

Eighteen  randomly  picked  points  are  treated  as  ob¬ 
servational  points  (Fig.  12).  The  reference  model  so¬ 
lution  at  the  18  points  is  the  observations,  O  =  (Ov  02 , 
....  018).  Using  the  optimization  method,  we  obtain  the 
temporally  varying  b , ( 1)  and  bit)  for  the  three  cases. 
For  case  1  (Fig.  13a),  both  b1  and  b,  fit  the  reference 
values  (Fig.  9c)  very  well  until  the  35th  day.  After  that 
day,  by  and  b2  change  rapidly  with  time  and  finally  blow 
up  at  the  45th  day.  For  case  2  (Fig.  13b),  both  bt  and 
b2  fit  the  reference  values  (Fig.  9c)  very  well.  For  case 
3  (Fig.  13c),  both  /?,  and  b,  are  very  close  to  the  ref¬ 
erence  values  (Fig.  9c)  with  small  perturbations.  Figure 
13  tells  us  that  smoothing  on  SB  is  a  key  issue  for  this 
method. 

We  integrate  the  POM  for  domain  A  with  the  com¬ 
puted  open  boundary  conditions  from  b,  and  b,  for  the 
three  cases.  The  30th  day’s  horizontal  velocity  fields 
(Fig.  14)  for  the  three  cases  all  agree  quite  well  with 
the  reference  model  results  (Fig.  11a),  and  the  90th  day’s 


horizontal  velocity  fields  (Fig.  15)  for  cases  2  and  3 
agree  quite  well  with  the  reference  model  results  (Fig. 
lib). 

e.  Error  estimation 

Similar  to  the  linear  case,  we  use  18  interior  obser¬ 
vations  O,  the  reference  boundary  vector  Bref  = 
bf),  the  open  boundary  vector  B  =  (/>,,  b2),  and  the 
interior  solution  S  =  (5),  S2,  .  .  . ,  Sis)  for  error  esti¬ 
mation. 

The  boundary  errors  E(B)  for  the  three  cases  are  shown 
in  Fig.  16.  When  there  is  no  smoothing  on  SB  (case  1), 
ElH'  keeps  very  small  values  (10  12-10  l  !)  for  the  first 
few  days.  After  the  eighth  day,  EiHl  increases  exponen¬ 
tially  with  time  and  reaches  the  order  of  1  at  the  45th 
day.  This  indicates  that  we  cannot  use  the  optimization 
method  without  smoothing  on  SB.  When  there  is 
smoothing  on  SB  and  no  noise  added  to  the  observations 
(case  2),  E(B)  has  larger  values  (~1(D2)  than  case  1  at 
the  beginning,  then  rapidly  decreases  with  time  during 
the  first  10  days  and  gradually  decreases  with  time  af¬ 
terward.  After  10  days  of  integration,  the  magnitude  of 
E(B)  is  on  the  order  of  10  4— 10  5.  When  there  is  smooth¬ 
ing  on  SB  and  noise  added  to  the  observations  (case  3), 
E(B)  fluctuates  around  10~25  (—3.162  X  1CD3)  with  the 
maximum  value  near  10  15  (—0.03)  and  the  minimum 
value  around  10  '.  Figure  16  indicates  that  smoothing 
on  SB  is  very  important  for  this  method. 

The  interior  errors  E(°>  for  the  three  cases  are  shown 
in  Fig.  17.  When  there  is  no  smoothing  on  SB  (case  1), 
E,a>  keeps  very  small  values  (10  l!)  for  the  first  few 
days  and  then  increases  exponentially  with  time  and 
reaches  the  order  of  1  at  the  45th  day.  When  there  is 
smoothing  on  SB  and  no  noise  added  to  the  observations 
(case  2),  Eto 1  has  larger  values  (~1(D2)  than  case  1  at 
the  beginning  then  rapidly  decreases  with  time  during 
first  10  days  and  gradually  decreases  with  time  after¬ 
ward.  After  10  days  of  integration,  the  magnitude  of 
E<0>  is  on  the  order  of  10  -10  \  When  there  is  smooth¬ 
ing  on  SB  and  noise  added  to  the  observations  (case  3), 
Em  fluctuates  around  1CD2-2  (—6.3  X  1CD3)  with  the 
maximum  value  near  1 0  2  and  the  minimum  value 
around  1CD24  (—3.9  X  1CD3).  Both  Figs.  16  and  17 
indicate  that  smoothing  on  SB  is  very  important  for 
determining  open  boundary  conditions. 

5.  Conclusions 

1)  The  proposed  optimization  method  provides  a  useful 
scheme  to  obtain  unknown  open  boundary  values 
from  known  interior  values.  Different  from  the  ad¬ 
joint  method,  this  scheme  can  be  easily  incorporated 
into  any  ocean  models. 

2)  Extremely  small  computational  errors  are  found  in 
applying  this  method  to  the  Csanady  shelf  model, 
which  proves  the  feasibility  of  using  this  optimiza¬ 
tion  method  for  linear  models. 
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3)  For  time-dependent  dynamical  models,  when  the 
temporally  varying  values  are  given  at  interior  ob¬ 
servation  points,  the  optimization  method  can  be 
used  for  each  time  step  to  obtain  the  unknown  open 
boundary  values  for  that  time  step. 

4)  For  a  primitive  equation  model  with  turbulent  mixing 
processes  (e.g.,  POM),  it  is  very  important  to  use 
smoothing  on  the  open  boundary  parameter  vector. 
If  smoothing  is  not  used,  POM  can  be  integrated 
only  within  a  certain  period  (45  days  in  our  case) 
and  will  blow  up  afterward.  If  smoothing  is  used, 
the  model  is  computationally  stable. 

5)  This  optimization  method  performs  well  even  when 
random  noises  are  added  to  the  observational  points 
(case  3).  This  indicates  that  we  can  use  real-time 
data  to  invert  for  the  unknown  open  boundary  values. 
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